﻿// gdal_app.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
// 使用rfef 文本文件生成 geotiff格式文件。波段从小到大表示高度层从低到高。

#include <iostream>
#include "gdal_priv.h"
#include <stdexcept>
using namespace std;
#include "RFEFAsciiReader.h"
#include "ajson5.h"
using namespace ArduinoJson;
#include "GIMReader.h"
#include <sstream>

int main(int argc, char* argv[])
{
    std::cout << "version v0.1" << endl ;
    std::cout << "version v0.2 add VTEC Volume Total Electron Concentration." << endl;
    std::cout << "version v0.3 add GIM" << endl;
    std::cout << "version v0.3.1 GIM 73x73" << endl;//为了适配在cesium显示，在纬度方向上下各增加一个像素凑成73x73，覆盖范围+/-90.0
    std::cout << "usage: elecden2tif mode /dir/filename.rfef /out/dir/output.tif" << endl;

    cout << "mode: rfef/vtec/gim" << endl;
    if (argc != 4) {
        return 11;
    }
    string mode = argv[1];
    string infilename = argv[2];
    string outfilename = argv[3];
    cout << "mode: " << mode << endl;
    std::cout << "in: " << infilename << endl;
    std::cout << "out: " << outfilename << endl;
    double scale = 10E7; // for rfef
    if (mode.compare("vtec") == 0) scale = 1;
    std::cout << "using scale: "<<scale << endl;
    GDALAllRegister();

    try {

        if (mode.compare("rfef") == 0 || mode.compare("vtec") == 0) {
            RFEFAsciiReader reader;
            bool asciiOk = reader.read(infilename);

            GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
            if (driver == NULL) {
                throw runtime_error("bad output file");
            }

            GDALDataset* outds = driver->Create(outfilename.c_str(), reader._dim1Size, reader._dim0Size, reader._dim2Size, GDT_Float32, NULL);
            if (outds == NULL) {
                throw runtime_error("failed to create output file.");
            }
            for (int ih = 0; ih < reader._dim2Size;++ih) {
                vector<float> heightData2D;
                heightData2D.resize(reader._dim0Size * reader._dim1Size);
                for (int iy = 0; iy < reader._dim0Size; ++iy) {
                    for (int ix = 0; ix < reader._dim1Size;++ix) {
                        double tval0 = reader.getValue3D(ix, iy, ih);
                        float tval1 = tval0 / scale;
                        //这里翻转，保证图像左上坐标对应纬度90
                        int iy_invert = reader._dim0Size - 1 - iy;
                        heightData2D[iy_invert * reader._dim1Size + ix] = tval1;
                    }
                }
                //write to disk
                outds->GetRasterBand(ih + 1)->RasterIO(GF_Write, 0, 0, reader._dim1Size, reader._dim0Size, heightData2D.data(),
                    reader._dim1Size, reader._dim0Size, GDT_Float32, 0, 0, 0);
            }
            double trans[6] = { -180, 360.0 / reader._dim1Size, 0.0 , 90.0 , 0.0, -180.0 / reader._dim0Size };
            //trans这里我使用常用的图像左上角对应 -180,90 地理坐标。像素数据我将Ascii的纬度顺序进行翻转了，以对应trans[5]的负值。
            outds->SetGeoTransform(trans);
            string wktStr = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,"
                "AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,"
                "AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],"
                "AUTHORITY[\"EPSG\",\"4326\"]]";
            outds->SetProjection(wktStr.c_str());

            GDALClose(outds);
        }
        else if (mode.compare("gim") == 0) {
            GIMReader reader;
            GIMReader::GIMData gimData;
            reader.read(infilename, gimData);
            cout << "gim " << gimData._valuesPerHour.size() << endl;
            GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
            if (driver == NULL) {
                throw runtime_error("bad output file");
            }
            for (int idt = 0; idt < gimData._datetimes.size(); ++idt) {
                std::stringstream ss;
                ss << outfilename << gimData._datetimes[idt] << ".tif";
                std::string outname1 = ss.str(); 
                GDALDataset* outds = driver->Create(outname1.c_str(),
                    gimData._nlon, gimData._nlat, 1 , GDT_Float32, NULL);
                if (outds == NULL) {
                    throw runtime_error("failed to create output file.");
                }

                //write to disk
                outds->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, gimData._nlon, gimData._nlat, gimData._valuesPerHour[idt].data(),
                    gimData._nlon, gimData._nlat, GDT_Float32, 0, 0, 0);
                
                double trans[6] = { gimData._lon1 , (gimData._lon2-gimData._lon1)/gimData._nlon, 0.0 , gimData._lat1 , 0.0,
                    (gimData._lat2-gimData._lat1) /gimData._nlat };
                outds->SetGeoTransform(trans);
                string wktStr = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,"
                    "AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,"
                    "AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],"
                    "AUTHORITY[\"EPSG\",\"4326\"]]";
                outds->SetProjection(wktStr.c_str());
                GDALClose(outds);
            }
            
        }
        
        //由于cesium读取tiff格式数据的代码已经验证通过了，所以后面应该不需要json数据了。
        if(false)
        {
            
        }


        cout << "done" << endl;
    }
    catch (exception& ex) {
        cout << "exception: " << ex.what() << endl;
        return 1;
    }
    return 0;
}

